Luis Alejandro Rodríguez Arenas Cod. 202321287¶

BLOQUE 1 Imports y configuración¶

In [1]:
import numpy as np
from numpy import linspace
import matplotlib.pyplot as plt
from sympy import symbols, sin, cos, pi, integrate, lambdify, simplify, latex, exp
from sympy import init_printing
from sympy import Piecewise
from matplotlib.animation import FuncAnimation

init_printing(use_latex='mathjax') 
In [2]:
x, t = symbols('x t', real=True)

BLOQUE 2 - Parámetros físicos del problema¶

ej 1 10-5¶

image.png

In [ ]:
L = 50
alpha = 2
T0 = 0
TL = 0

f = 20

print("Función inicial f(x) =")
display(f)

ej 1 10-6¶

image.png

In [3]:
L = 30
alpha = 2
T0 = 20
TL = 50

f = 60-2*x

print("Función inicial f(x) =")
display(f)
Función inicial f(x) =
$\displaystyle 60 - 2 x$

ej 1 10-7¶

image.png

In [16]:
L = 30
a = 4
T0 = 0
TL = 0 

f = Piecewise(
    (x/10,  (0 <= x) & (x <= 10)),
    ((30-x)/20, (10 < x) & (x <= 30)),
    (0, x > 30),
)

print("Función inicial f(x) por tramos =")
display(f)
Función inicial f(x) por tramos =
$\displaystyle \begin{cases} \frac{x}{10} & \text{for}\: x \geq 0 \wedge x \leq 10 \\\frac{3}{2} - \frac{x}{20} & \text{for}\: x \leq 30 \wedge x > 10 \\0 & \text{for}\: x > 30 \end{cases}$

Ejemplo para molestar a FTCS¶

In [ ]:
L = 30   
alpha = 2
T0 = 0
TL = 0

f = sin(15*pi*x/L) + sin(35*pi*x/L)

Separación en parte estacionaria + transitoria¶

In [17]:
# Parte estacionaria v(x): resolviendo v''(x)=0, v(0)=T0, v(L)=TL
v = T0 + (TL - T0) * (x/L)

print("Solución estacionaria v(x) =")
display(v)

# Condición inicial para v(x,0): v0(x) = f(x) – s(x)
v0 = simplify(f - v)

print("\nv(x,0) = f(x) - s(x) =")
display(v0)
Solución estacionaria v(x) =
$\displaystyle 0$
v(x,0) = f(x) - s(x) =
$\displaystyle \begin{cases} \frac{x}{10} & \text{for}\: x \geq 0 \wedge x \leq 10 \\\frac{3}{2} - \frac{x}{20} & \text{for}\: x \leq 30 \wedge x > 10 \\0 & \text{for}\: x > 30 \\\text{NaN} & \text{otherwise} \end{cases}$

BLOQUE 4 - Serie de Fourier¶

image-3.png image.png

In [18]:
# Fórmula simbólica del coeficiente c_n
n = symbols('n', integer=True, positive=True)
cn_formula = (2/L) * integrate(v0 * sin(n*pi*x/L), (x, 0, L))

print("Coeficiente general c_n =")
display(simplify(cn_formula))
Coeficiente general c_n =
$\displaystyle \frac{9.0 \sin{\left(\frac{\pi n}{3} \right)}}{\pi^{2} n^{2}}$
In [6]:
# Para calor
print("u(x,T) =")
display(cn_formula * exp(-(alpha*n*pi/L)**2 * t) * sin(n*pi*x/L))
u(x,T) =
$\displaystyle \left(\frac{100.0 \left(-1\right)^{n}}{\pi n} + \frac{80.0}{\pi n}\right) e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}$
In [19]:
# Para onda
print("u(x,t) =")
display(cn_formula * cos(-(alpha*n*pi/L) * t) * sin(n*pi*x/L))
u(x,t) =
$\displaystyle \frac{9.0 \sin{\left(\frac{\pi n}{3} \right)} \sin{\left(\frac{\pi n x}{30} \right)} \cos{\left(\frac{\pi n t}{15} \right)}}{\pi^{2} n^{2}}$

Cálculo de coeficientes numéricos (Calor)¶

In [7]:
N = 100
cn = []

for k in range(1, N+1):
    val = float(cn_formula.subs(n, k))
    cn.append(val)

print("\nPrimeros coeficientes c_n (numéricos):")
for i, c in enumerate(cn[:8], 1):
    print(f"c{i} = {c:.6f}")

u_series_sym = 0
for k in range(1, N+1):
    ck = cn_formula.subs(n, k)
    term = ck * exp(-(alpha*n*pi/L)**2 * t) * sin(n*pi*x/L)
    u_series_sym += term

print("\nSolución v(x,t) truncada simbólica (antes de sumar v(x)):")
display(simplify(u_series_sym))

u_full_sym = v + u_series_sym

print("\nSolución completa u(x,t)")
display(simplify(u_full_sym))
Primeros coeficientes c_n (numéricos):
c1 = -6.366198
c2 = 28.647890
c3 = -2.122066
c4 = 14.323945
c5 = -1.273240
c6 = 9.549297
c7 = -0.909457
c8 = 7.161972

Solución v(x,t) truncada simbólica (antes de sumar v(x)):
$\displaystyle \frac{346.17298348015 e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}}{\pi}$
Solución completa u(x,t)
$\displaystyle x + 20 + \frac{346.17298348015 e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}}{\pi}$
Cálculo de coeficientes numéricos (Onda)¶
In [20]:
# Coeficiente correcto para ONDA
A_n_formula = (2/L) * integrate(f * sin(n*pi*x/L), (x, 0, L))

N = 100
A = []

for k in range(1, N+1):
    A.append(float(A_n_formula.subs(n, k)))

print("Coeficientes A_n de la onda:")
for i, c in enumerate(A[:8], 1):
    print(f"A{i} = {c:.6f}")

def u_series_onda(x_array, t_value, A, N):
    total = np.zeros_like(x_array)
    for k in range(1, N+1):
        total += A[k-1] * np.sin(k*np.pi*x_array/L) * np.cos((k*np.pi*a/L)*t_value)
    return total

def u_sol_onda(x_array, t_value):
    return u_series_onda(x_array, t_value, A, N)
Coeficientes A_n de la onda:
A1 = 0.789720
A2 = 0.197430
A3 = 0.000000
A4 = -0.049358
A5 = -0.031589
A6 = 0.000000
A7 = 0.016117
A8 = 0.012339

BLOQUE 5 - Construcción de u(x,t)¶

Calor¶

In [8]:
x_vals = np.linspace(0, L, 300)

def v_series_calor(x_array, t_value, cn, N):
    total = np.zeros_like(x_array)
    for k in range(1, N+1):
        total += cn[k-1] * np.sin(k*np.pi*x_array/L) * np.exp(-alpha*(k*np.pi/L)**2 * t_value)
    return total

v_func = lambdify(x, v, 'numpy')

def u_sol_calor(x_array, t_value):
    return v_func(x_array) + v_series_calor(x_array, t_value, cn, N)

u_sol_calor(x_vals, 0)
Out[8]:
array([20.        , 45.21579184, 61.84221647, 66.55056011, 62.68915318,
       57.07310592, 54.92216581, 56.86623359, 59.89584709, 60.82157044,
       59.04915746, 56.5541714 , 55.61603767, 56.6844463 , 58.25590197,
       58.56542695, 57.28372131, 55.6454408 , 55.0841114 , 55.83566643,
       56.83905447, 56.89287802, 55.83423964, 54.60252645, 54.22685952,
       54.81180989, 55.50833404, 55.42075765, 54.49207593, 53.50127626,
       53.23719968, 53.71732452, 54.21885641, 54.04075526, 53.19892705,
       52.36984744, 52.18083512, 52.58739073, 52.95208813, 52.71062262,
       51.93227901, 51.22093087, 51.08620101, 51.43709715, 51.69904853,
       51.41050338, 50.68163564, 50.06106032, 49.96759107, 50.27397334,
       50.45487023, 50.12984184, 49.44145415, 48.89393664, 48.83297236,
       49.10219033, 49.21668866, 48.86251043, 48.20853564, 47.72181226,
       47.68713577, 47.92424236, 47.9827078 , 47.60470862, 46.98090689,
       46.54613737, 46.53313528, 46.7417106 , 46.75174329, 46.35395393,
       45.75728636, 45.36788925, 45.37300994, 45.55564437, 45.52298103,
       45.10855567, 44.53680701, 44.18775212, 44.20817365, 44.36676601,
       44.29584122, 43.86732205, 43.31886259, 43.00622109, 43.03963833,
       43.17558768, 43.06989791, 42.62938806, 42.10301784, 41.82366521,
       41.86814803, 41.98248127, 41.84482917, 41.39410991, 40.88895354,
       40.64036724, 40.69426294, 40.78772192, 40.62038499, 40.1609976 ,
       39.67643169, 39.45654966, 39.5184138 , 39.59151611, 39.39636605,
       38.92967062, 38.46527275, 38.27239229, 38.34093835, 38.39402021,
       38.172609  , 37.69982784, 37.25534039, 37.08804454, 37.1621063 ,
       37.19535311, 36.94897619, 36.47122657, 36.046531  , 35.90363422,
       35.98213698, 35.99560492, 35.725348  , 35.24366769, 34.8387665 ,
       34.71927417, 34.80121207, 34.79484297, 34.5016171 , 34.01698465,
       33.63198922, 33.53506735, 33.61948496, 33.59311609, 33.27768389,
       32.79103541, 32.42615846, 32.35111096, 32.43708793, 32.3904575 ,
       32.05345266, 31.5656961 , 31.22124803, 31.16750005, 31.25413769,
       31.18688673, 30.82882821, 30.34085603, 30.01724476, 29.98433075,
       30.07073989, 29.98241075, 29.60371268, 29.11641364, 28.81414773,
       28.80170349, 28.88699291, 28.7770245 , 28.37800228, 27.8922729 ,
       27.61196805, 27.6197263 , 27.70299124, 27.57071073, 27.15158377,
       26.66834008, 26.41072929, 26.43851859, 26.51882863, 26.36343927,
       25.92433049, 25.44452061, 25.21046855, 25.25821546, 25.3346013 ,
       25.15516562, 24.69609759, 24.22071575, 24.01123828, 24.07897303,
       24.15041139, 23.94582874, 23.46671602, 22.99681899, 22.81310894,
       22.90097513, 22.96637089, 22.7353478 , 22.23598499, 21.77271169,
       21.61617286, 21.72444203, 21.78260639, 21.52361756, 21.00366187,
       20.54825777, 20.42054978, 20.54964203, 20.59926517, 20.31050192,
       19.7694488 , 19.32329682, 19.2263945 , 19.37690719, 19.4165231 ,
       19.09582465, 18.53297408, 18.09763481, 18.03390789, 18.20665526,
       18.23459547, 17.87935626, 17.29376619, 16.87103144, 16.84335279,
       17.03942085, 17.05375213, 16.66079485, 16.05121621, 15.64318192,
       15.65507726, 15.87590085, 15.87433931, 15.43973752, 14.80452246,
       14.41369056, 14.46954978, 14.71702246, 14.69681203, 14.21563679,
       13.55260642, 13.1820308 , 13.28741371, 13.56404803, 13.5217837 ,
       12.98773191, 12.29398071, 11.94748305, 12.10957474, 12.41874273,
       12.35010478, 11.75493651, 11.02653404, 10.70903413, 10.93734691,
       11.28365344, 11.1829927 , 10.51564645,  9.74716512,  9.46520749,
        9.77270872, 10.1625959 , 10.02225659,  9.2673933 ,  8.45112487,
        8.21376146,  8.61877974,  9.06155727,  8.87070887,  8.00617681,
        7.13075527,  6.95111824,  7.48077522,  7.99049502,  7.7329729 ,
        6.7250672 ,  5.77286293,  5.67119607,  6.36810831,  6.96727509,
        6.61721362,  5.41094316,  4.35262857,  4.36277122,  5.29964075,
        6.0274404 ,  5.5392927 ,  4.03566136,  2.81724289,  3.00267863,
        4.31934875,  5.25311901,  4.53440507,  2.52634234,  1.03131863,
        1.5348144 ,  3.55779259,  4.88545709,  3.69735932,  0.62219005,
       -1.48410087, -0.21203353,  3.6298674 ,  6.05004443,  3.34766987,
       -3.62505196, -8.34057928, -2.34202258, 18.49211875, 50.        ])

Onda¶

In [21]:
x_vals = np.linspace(0, L, 300)

def u_series_onda(x_array, t_value, A, N):
    total = np.zeros_like(x_array)
    for k in range(1, N+1):
        total += A[k-1] * np.sin(k*np.pi*x_array/L) * np.cos((k*np.pi*a/L)*t_value)
    return total

def u_sol_onda(x_array, t_value):
    return u_series_onda(x_array, t_value, A, N)

u_sol_onda(x_vals, 0)
Out[21]:
array([0.00000000e+00, 1.00343941e-02, 2.00690687e-02, 3.01027554e-02,
       4.01339423e-02, 5.01636934e-02, 6.01958188e-02, 7.02329527e-02,
       8.02725683e-02, 9.03082994e-02, 1.00336495e-01, 1.10361270e-01,
       1.20391570e-01, 1.30431962e-01, 1.40476607e-01, 1.50513945e-01,
       1.60538600e-01, 1.70558272e-01, 1.80587183e-01, 1.90631421e-01,
       2.00681262e-01, 2.10719770e-01, 2.20740249e-01, 2.30754609e-01,
       2.40782575e-01, 2.50831347e-01, 2.60886641e-01, 2.70925863e-01,
       2.80941413e-01, 2.90950150e-01, 3.00977646e-01, 3.11031783e-01,
       3.21092900e-01, 3.31132335e-01, 3.41142033e-01, 3.51144712e-01,
       3.61172268e-01, 3.71232804e-01, 3.81300260e-01, 3.91339335e-01,
       4.01342015e-01, 4.11338027e-01, 4.21366273e-01, 4.31434532e-01,
       4.41509041e-01, 4.51547061e-01, 4.61541205e-01, 4.71529706e-01,
       4.81559421e-01, 4.91637162e-01, 5.01719725e-01, 5.11755801e-01,
       5.21739354e-01, 5.31719152e-01, 5.41751362e-01, 5.51841011e-01,
       5.61933058e-01, 5.71965986e-01, 5.81936050e-01, 5.91905415e-01,
       6.01941558e-01, 6.12046621e-01, 6.22150267e-01, 6.32178300e-01,
       6.42130570e-01, 6.52086887e-01, 6.62129129e-01, 6.72254978e-01,
       6.82373501e-01, 6.92393894e-01, 7.02321538e-01, 7.12260628e-01,
       7.22312536e-01, 7.32468056e-01, 7.42606892e-01, 7.52614871e-01,
       7.62506024e-01, 7.72420643e-01, 7.82488836e-01, 7.92690392e-01,
       8.02859483e-01, 8.12845460e-01, 8.22676624e-01, 8.32552620e-01,
       8.42651779e-01, 8.52934910e-01, 8.63155356e-01, 8.73095134e-01,
       8.82808667e-01, 8.92612615e-01, 9.02787010e-01, 9.13255031e-01,
       9.23584207e-01, 9.33382600e-01, 9.42761832e-01, 9.52384821e-01,
       9.62908697e-01, 9.74172342e-01, 9.84799512e-01, 9.92665184e-01,
       9.96028152e-01, 9.94574765e-01, 9.89622574e-01, 9.83327655e-01,
       9.77464303e-01, 9.72626103e-01, 9.68324356e-01, 9.63758414e-01,
       9.58575439e-01, 9.53061104e-01, 9.47744872e-01, 9.42875010e-01,
       9.38241973e-01, 9.33445966e-01, 9.28297652e-01, 9.22967190e-01,
       9.17777975e-01, 9.12881337e-01, 9.08131037e-01, 9.03254207e-01,
       8.98126194e-01, 8.92877234e-01, 8.87744197e-01, 8.82831893e-01,
       8.78021509e-01, 8.73103193e-01, 8.67987539e-01, 8.62784821e-01,
       8.57683811e-01, 8.52761552e-01, 8.47914387e-01, 8.42970492e-01,
       8.37863006e-01, 8.32690319e-01, 8.27610293e-01, 8.22681225e-01,
       8.17809012e-01, 8.12847523e-01, 8.07745794e-01, 8.02594371e-01,
       7.97529366e-01, 7.92595370e-01, 7.87704831e-01, 7.82730322e-01,
       7.77632859e-01, 7.72497435e-01, 7.67443863e-01, 7.62506136e-01,
       7.57601473e-01, 7.52616811e-01, 7.47522639e-01, 7.42399810e-01,
       7.37355343e-01, 7.32414686e-01, 7.27498694e-01, 7.22505797e-01,
       7.17414254e-01, 7.12301697e-01, 7.07264732e-01, 7.02321701e-01,
       6.97396329e-01, 6.92396549e-01, 6.87307168e-01, 6.82203230e-01,
       6.77172618e-01, 6.72227605e-01, 6.67294265e-01, 6.62288592e-01,
       6.57201039e-01, 6.52104505e-01, 6.47079388e-01, 6.42132678e-01,
       6.37192424e-01, 6.32181605e-01, 6.27095636e-01, 6.22005589e-01,
       6.16985310e-01, 6.12037112e-01, 6.07090747e-01, 6.02075363e-01,
       5.96990800e-01, 5.91906532e-01, 5.86890573e-01, 5.81941041e-01,
       5.76989191e-01, 5.71969706e-01, 5.66886416e-01, 5.61807371e-01,
       5.56795317e-01, 5.51844563e-01, 5.46887724e-01, 5.41864513e-01,
       5.36782400e-01, 5.31708134e-01, 5.26699644e-01, 5.21747751e-01,
       5.16786321e-01, 5.11759695e-01, 5.06678689e-01, 5.01608844e-01,
       4.96603633e-01, 4.91650660e-01, 4.86684963e-01, 4.81655182e-01,
       4.76575233e-01, 4.71509517e-01, 4.66507345e-01, 4.61553333e-01,
       4.56583634e-01, 4.51550920e-01, 4.46471996e-01, 4.41410167e-01,
       4.36410830e-01, 4.31455804e-01, 4.26482322e-01, 4.21446866e-01,
       4.16368948e-01, 4.11310806e-01, 4.06314124e-01, 4.01358099e-01,
       3.96381017e-01, 3.91342985e-01, 3.86266064e-01, 3.81211442e-01,
       3.76217260e-01, 3.71260239e-01, 3.66279710e-01, 3.61239250e-01,
       3.56163325e-01, 3.51112084e-01, 3.46120263e-01, 3.41162242e-01,
       3.36178394e-01, 3.31135637e-01, 3.26060716e-01, 3.21012737e-01,
       3.16023153e-01, 3.11064123e-01, 3.06077065e-01, 3.01032127e-01,
       2.95958223e-01, 2.90913407e-01, 2.85925949e-01, 2.80965894e-01,
       2.75975715e-01, 2.70928704e-01, 2.65855835e-01, 2.60814098e-01,
       2.55828665e-01, 2.50867563e-01, 2.45874342e-01, 2.40825354e-01,
       2.35753543e-01, 2.30714816e-01, 2.25731314e-01, 2.20769141e-01,
       2.15772942e-01, 2.10722065e-01, 2.05651339e-01, 2.00615563e-01,
       1.95633908e-01, 1.90670635e-01, 1.85671509e-01, 1.80618826e-01,
       1.75549217e-01, 1.70516343e-01, 1.65536456e-01, 1.60572049e-01,
       1.55570043e-01, 1.50515629e-01, 1.45447171e-01, 1.40417159e-01,
       1.35438967e-01, 1.30473391e-01, 1.25468538e-01, 1.20412465e-01,
       1.15345194e-01, 1.10318015e-01, 1.05341448e-01, 1.00374665e-01,
       9.53669936e-02, 9.03093269e-02, 8.52432843e-02, 8.02189119e-02,
       7.52439069e-02, 7.02758747e-02, 6.52654056e-02, 6.02062077e-02,
       5.51414362e-02, 5.01198537e-02, 4.51463503e-02, 4.01770244e-02,
       3.51637715e-02, 3.01031008e-02, 2.50396463e-02, 2.00208429e-02,
       1.50487844e-02, 1.00781174e-02, 5.06208853e-03, 2.67536267e-16])

BLOQUE 6 - Gráficas¶

Calor¶

In [9]:
x_vals = np.linspace(0, L, 300)
t_list = [0, 1, 10, 100, 500]

plt.figure(figsize=(8,5))

for ti in t_list:
    ui = u_sol_calor(x_vals, ti)
    plt.plot(x_vals, ui, label=f"t={ti}s")

plt.title("Solución analítica por Serie de Fourier")
plt.xlabel("x")
plt.ylabel("Temperatura u(x,t)")
plt.legend()
plt.grid(True)
plt.show()
No description has been provided for this image

Onda¶

In [ ]:
x_vals = np.linspace(0, L, 300)
t_list = [0, 1, 10, 100, 500]

plt.figure(figsize=(8,5))

for ti in t_list:
    ui = u_sol_onda(x_vals, ti)
    plt.plot(x_vals, ui, label=f"t={ti}s")

plt.title("Solución analítica por Serie de Fourier")
plt.xlabel("x")
plt.ylabel("Temperatura u(x,t)")
plt.legend()
plt.grid(True)
plt.show()
No description has been provided for this image

BLOQUE 7 - Animación¶

Calor¶

In [10]:
fig, ax = plt.subplots(figsize=(8,5))

line, = ax.plot([], [], lw=2)

ax.set_xlim(0, L)
ax.set_ylim(np.min(u_sol_calor(x_vals, 500)), np.max(u_sol_calor(x_vals, 0)))
ax.grid(True)

t_anim = np.linspace(0, 200, 300)

def init():
    line.set_data([], [])
    return line,

def update(frame):
    t_val = t_anim[frame]
    u_val = u_sol_calor(x_vals, t_val)
    line.set_data(x_vals, u_val)
    ax.set_title(f"t = {t_val:.1f} s")
    return line,

anim = FuncAnimation(
    fig,
    update,
    frames=len(t_anim),
    init_func=init,
    interval=100,
    blit=False
)


from IPython.display import HTML
HTML(anim.to_jshtml())
Out[10]:
No description has been provided for this image
No description has been provided for this image

Onda¶

In [23]:
fig, ax = plt.subplots(figsize=(8,5))

line, = ax.plot([], [], lw=2)

modelo = "onda"

if modelo == "calor":
    u_test_min = u_sol_calor(x_vals, 500)
    u_test_max = u_sol_calor(x_vals, 0)
    ax.set_ylim(np.min(u_test_min), np.max(u_test_max))

elif modelo == "onda":
    u_test = u_sol_onda(x_vals, 0)
    ax.set_ylim(np.min(u_test)*1.1, np.max(u_test)*1.1)

ax.set_xlim(0, L)
ax.grid(True)

if modelo == "calor":
    t_anim = np.linspace(0, 500, 300)
else:
    t_anim = np.linspace(0, 10, 300)

def init():
    line.set_data([], [])
    return line,

def update(frame):
    t_val = t_anim[frame]
    
    if modelo == "calor":
        u_val = u_sol_calor(x_vals, t_val)
    else:
        u_val = u_sol_onda(x_vals, t_val)
    
    line.set_data(x_vals, u_val)
    ax.set_title(f"t = {t_val:.2f} s")
    return line,

anim = FuncAnimation(
    fig,
    update,
    frames=len(t_anim),
    init_func=init,
    interval=50,
    blit=False
)

from IPython.display import HTML
HTML(anim.to_jshtml())
Out[23]:
No description has been provided for this image
No description has been provided for this image

BLOQUE 8 - Superficie 3D u(x,t)¶

Calor¶

In [11]:
from mpl_toolkits.mplot3d import Axes3D

# Mallado de espacio y tiempo
x_vals = np.linspace(0, L, 200)
t_vals = np.linspace(0, 500, 1000)

X, T = np.meshgrid(x_vals, t_vals)

# Matriz U(x,t)
U = np.zeros_like(X)

for i in range(len(t_vals)):
    U[i, :] = u_sol_calor(x_vals, t_vals[i])

# Figura 3D
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
# Rotar la vista 3D
ax.view_init(elev=30, azim=45)

surf = ax.plot_surface(X, T, U, cmap='inferno', edgecolor='none')

ax.set_title("Evolución térmica u(x,t) — Ecuación de Calor")
ax.set_xlabel("Posición x")
ax.set_ylabel("Tiempo t")
ax.set_zlabel("Temperatura u")

fig.colorbar(surf, shrink=0.5, aspect=8)
plt.show()
No description has been provided for this image

Onda¶

In [24]:
from mpl_toolkits.mplot3d import Axes3D

# Mallado de espacio y tiempo
x_vals = np.linspace(0, L, 200)
t_vals = np.linspace(0, 80, 200)  # la onda oscila rápido → tiempos pequeños

X, T = np.meshgrid(x_vals, t_vals)

# Matriz U(x,t)
U = np.zeros_like(X)

for i in range(len(t_vals)):
    U[i, :] = u_sol_onda(x_vals, t_vals[i])

# Figura 3D
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
# Rotar la vista 3D
ax.view_init(elev=30, azim=45)
surf = ax.plot_surface(X, T, U, cmap='viridis', edgecolor='none')

ax.set_title("Evolución u(x,t) — Ecuación de Onda")
ax.set_xlabel("Posición x")
ax.set_ylabel("Tiempo t")
ax.set_zlabel("Desplazamiento u")

fig.colorbar(surf, shrink=0.5, aspect=8)
plt.show()
No description has been provided for this image

Bloque 9 - Solución con FTCS y comparaciones¶

Calor¶

In [14]:
def explicit_fd(f_init, L, a, T0, TL, nx, t_max):

    dx = L/(nx-1)
    
    # dt controlado por estabilidad CFL
    dt = 0.45 * dx**2 / a     
    nt = int(t_max / dt) + 1

    x_grid = np.linspace(0, L, nx)
    times = np.linspace(0, t_max, nt)

    u = f_init(x_grid)
    u[0], u[-1] = T0, TL
    
    hist = np.zeros((nt, nx))
    hist[0] = u.copy()

    r = a * dt / dx**2

    for k in range(1, nt):
        un = u.copy()
        u[1:-1] = un[1:-1] + r*(un[2:] - 2*un[1:-1] + un[:-2])
        u[0], u[-1] = T0, TL
        hist[k] = u.copy()

    return x_grid, times, hist

nx_num = 10
t_max  = 80.0

print(f"Ejecutando FTCS con nx = {nx_num}")

f_init_fn = lambda X: u_sol_calor(X, 0)

x_num, t_num, U_num_hist = explicit_fd(
    f_init_fn, L, alpha, T0, TL,
    nx=nx_num,
    t_max=t_max
)

U_analytic_hist = np.zeros_like(U_num_hist)
for i, tt in enumerate(t_num):
    U_analytic_hist[i] = u_sol_calor(x_num, tt)

t_compare = [0, 1, 10, 40]

plt.figure(figsize=(9,6))
for tt in t_compare:
    idx = np.argmin(np.abs(t_num - tt))
    plt.plot(x_num, U_analytic_hist[idx], lw=2, label=f"Analítica t={t_num[idx]:.2f}")
    plt.plot(x_num, U_num_hist[idx], '--', lw=1.5, label=f"FTCS t={t_num[idx]:.2f}")
plt.title(f"Comparación Analítica vs FTCS — nx={nx_num}")
plt.xlabel("x"); plt.ylabel("u(x,t)")
plt.grid(True); plt.legend()
plt.show()

x_surf, t_surf = x_num, t_num
X_surf, T_surf = np.meshgrid(x_surf, t_surf)

U_analytic_surf = U_analytic_hist
U_num_surf = U_num_hist

fig = plt.figure(figsize=(18,5))

# analítica
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X_surf, T_surf, U_analytic_surf, cmap='inferno')
ax1.set_title("Analítica u(x,t)")
ax1.view_init(30,45)

# numérica
ax2 = fig.add_subplot(132, projection='3d')
ax2.plot_surface(X_surf, T_surf, U_num_surf, cmap='inferno')
ax2.set_title("FTCS u_num(x,t)")
ax2.view_init(30,45)

# error
ax3 = fig.add_subplot(133, projection='3d')
ax3.plot_surface(X_surf, T_surf, U_num_surf - U_analytic_surf, cmap='seismic')
ax3.set_title("Error u_num - u_analytic")
ax3.view_init(30,45)

plt.tight_layout()
plt.show()
Ejecutando FTCS con nx = 10
No description has been provided for this image
No description has been provided for this image

Onda¶

In [26]:
nx_num = 10
t_max  = 10.0

f_init_fn = lambda X: u_sol_onda(X, 0)


x_num, t_num, U_num_hist = explicit_fd(
    f_init_fn, L, a, T0, TL,
    nx=nx_num,
    t_max=t_max
)

U_analytic_hist = np.zeros_like(U_num_hist)
for i, tt in enumerate(t_num):
    U_analytic_hist[i, :] = u_sol_onda(x_num, tt)


t_compare = [0.0, 1.0, 2.0, 4.0]

plt.figure(figsize=(9,6))
for tt in t_compare:
    idx = np.argmin(np.abs(t_num - tt))
    plt.plot(x_num, U_analytic_hist[idx], '-', lw=2, label=f'Analítica t={t_num[idx]:.2f}')
    plt.plot(x_num, U_num_hist[idx], '--', lw=1.5, label=f'FTCS t={t_num[idx]:.2f}')
plt.xlabel('x'); plt.ylabel('u(x,t)')
plt.title(f'Comparación Analítica vs FTCS (onda) — nx={nx_num}')
plt.legend(); plt.grid(True)
plt.show()

x_surf, t_surf, U_num_surf = explicit_fd(f_init_fn, L, a, T0, TL, nx=nx_num, t_max=t_max)

U_analytic_surf = np.zeros_like(U_num_surf)
for i, tt in enumerate(t_surf):
    U_analytic_surf[i, :] = u_sol_onda(x_surf, tt)

X_surf, T_surf = np.meshgrid(x_surf, t_surf)

fig = plt.figure(figsize=(18,5))

ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X_surf, T_surf, U_analytic_surf, cmap='viridis', edgecolor='none')
ax1.set_title("Analítica onda u(x,t)")
ax1.set_xlabel("x"); ax1.set_ylabel("t"); ax1.set_zlabel("u")
ax1.view_init(elev=30, azim=45)

ax2 = fig.add_subplot(132, projection='3d')
ax2.plot_surface(X_surf, T_surf, U_num_surf, cmap='viridis', edgecolor='none')
ax2.set_title("Numérica (FTCS) — NO exacta")
ax2.set_xlabel("x"); ax2.set_ylabel("t"); ax2.set_zlabel("u")
ax2.view_init(elev=30, azim=45)

ax3 = fig.add_subplot(133, projection='3d')
ax3.plot_surface(X_surf, T_surf, U_num_surf - U_analytic_surf, cmap='seismic', edgecolor='none')
ax3.set_title("Error FTCS - Analítica")
ax3.set_xlabel("x"); ax3.set_ylabel("t"); ax3.set_zlabel("error")
ax3.view_init(elev=30, azim=45)

plt.tight_layout()
plt.show()
No description has been provided for this image
No description has been provided for this image

Bloque 10 - Análisis de error y convergencia¶

In [15]:
Error = U_num_hist - U_analytic_hist

L2 = np.sqrt(np.sum(Error**2, axis=1) / Error.shape[1])
Linf = np.max(np.abs(Error), axis=1)

plt.figure(figsize=(8,4))
plt.plot(t_num, L2, label='L2 (RMS) del error')
plt.plot(t_num, Linf, label='L∞ del error')
plt.xlabel('Tiempo t'); plt.ylabel('Error')
plt.title('Normas del error en el tiempo')
plt.legend(); plt.grid(True)
plt.show()

plt.figure(figsize=(10,5))
plt.contourf(x_num, t_num, Error, levels=60, cmap='seismic')
plt.colorbar(label='u_num - u_analytic')
plt.xlabel('x'); plt.ylabel('t')
plt.title('Mapa de calor del error E(x,t)')
plt.show()

t_final = min(40.0, t_num[-1])
nx_list = [50, 100, 200, 400]
errors_final = []

for nx_test in nx_list:
    x_test, t_test, U_test = explicit_fd(f_init_fn, L, alpha, T0, TL, nx=nx_test, t_max=t_final)
    U_num_final = U_test[-1, :]
    U_analytic_final = u_sol_calor(x_test, t_final)
    err_L2 = np.sqrt(np.sum((U_num_final - U_analytic_final)**2) / U_num_final.size)
    errors_final.append(err_L2)

print("nx \t L2 error @ t_final")
for nxv, ev in zip(nx_list, errors_final):
    print(f"{nxv}\t {ev:.6e}")
No description has been provided for this image
No description has been provided for this image
nx 	 L2 error @ t_final
50	 4.219587e-03
100	 1.229874e-03
200	 1.882337e-04
400	 6.767316e-05

image.png image-2.png